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Abstract 

We present a transfer matrix method which is particularly useful for solving 
some classes of sandpile models. The method is then used to solve the de- 
terministic nonabelian sandpile models for N=2 and N=3. The possibility of 
generalization to arbitrary N is discussed briefly. 
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I. INTRODUCTION 



Sandpile model is the simplest model of self-organized criticality(SOC). This kind of 
models were first discussed by Bak, Tang, and Wiensenfeld They used numerical method 
to study some of the sandpile models and observed that the models automatically evolve into 
a self-organized critical state when they possess 1// spectra in both spatial and temporal 
distributions of certain physical quantities. They suspected that SOC maybe an universal 
characteristic underlying the the nonlinear dispersive systems such as earth quakes, forest 
fire, turbulence, etc, which are prevalent in nature. 

The simplest sandpile model is a cellular automata on an one dimensional lattice with 
a height number hi assigned to each site. There are two basic operations of the model - 
dropping and toppling. Dropping means that one sand is added at some site of the lattice, 
i.e. Di : hi ^ hi^i. Toppling occurs when a slope (defined as the difference in height 
between adjacent sites) exceeds some critical value. If toppling occurs at one site, some 
sands at the site will be moved to the other sites which may trigger further topplings. 

Though the rules of evolution of sandpile model are typically very simple, it is hard to 
solve them analytically when the degrees of freedom become very large. Most researcher 
handle the models by numerical simulation [0 ||^ If the rules are such that the evolution 
of the system is independent of the order of the droppings, then the model is called abelian. 
For a large class of abelian models some exact results have been obtained by Dhar et al. |^ 
The non-abelian ones are harder to solve and there exists little exact result. In a previous 
work 1^, one of the authors investigated a class of non-abelian sandpile models and was 
able to solve the model in the deterministic case when the sand is dropped at a fixed site. 
In this note, we wish to present a new method of solving this class of sandpile models. We 
consider the one dimensional case and label the sites from left to right as 1 to L. The sand 
is dropped only at the site 1. If the slope at a site exceeds a given number A^, then the sand 
will topple to the right. Let the slope cTi = hi ~ /ij+i, then the toppling rule is "if ai > N, 
then aj_i + N , ^ ai — {N + 1) , and cXj+Ar — > cTj+at + 1". The rule should be 
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modified when toppling occurs near the boundary. The condition at the left boundary is 
trivial. When sands reach beyond the right boundary they drop out from the system, i.e. 
we keep hi = for i > L. A state in which all cr, < is called a stable state. Toppling 
stops when a stable state is reached. Each dropping and subsequent toppling processes will 
result in the transition from one stable state to another. Since there are only a finite number 
of different states in the system, after dropping enough sand at site 1, the system will step 
into a cycle called the limit cycle. For the system that we consider here there is only one 
limit cycle in the problem. The number of different states in the limit cycle is A^^. In the 
following, all the states we will refer to are the stable states. 

The method to be introduced here is particularly useful for those models in which the 
structure of the limit cycles has been work out. Once the cycle structure is known, its 
information can be succinctly summarized in a matrix which we call the transfer matrix 
by analogy with the similar matrix in statistical mechanics. For the deterministic model 
defined above, the structure of the limit cycle has been worked out in ref. Therefore we 
should use it as our main example, the method may be applicable to a much wider classes 
of models. 

II. THE DEFINITION OF TRANSFER MATRIX 

A state is in the limit circle if and only if 

1. There exists at least one site i for any consecutive N sites such that CTj = A^. 

2. There exists a site i satisfying L — < i < L such that (Tj = N. 

Since the first condition is for consecutive sites, there is no constraint in consecutive 
A^ — 1 sites. Choose the complete set of states for consecutive A^ — 1 sites as a basis. For 
example in the N = 2 model, the states in the basis is |0), |2). For A^ = 3 model, the 
basis is of the form |33), |a3), |3a), \ab), where a, b belongs to {0, 1,2}. For arbitrary A^, 
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{Vi\T\Vi+i) 



there are (A^ + 1)^^"^^ states in the basis. Define the transfer matrix as a mapping between 
\r]i) = |(Ti(Ti+i • ■ ■ cri+Ar_2) and \r]i+i) = |ai+icri+2 ■ • ■ ai+N-i) such that 

1 if the sequence cTjCrj + 1 ■ ■ ■ (Tj+Ar_i satisfies constrain 1 ^ ^ 

otherwise 

For example, in the N = 2 model, if cxj+i = 2, the cxj could be 0, 1, 2, so (0|T|2) = 
(l|r|2) = (2|r|2) = 1. liai = 2, then a^+i could be 0, 1, ,2, so that (2|T|0) = (2|r|l) = 1. 
For other cases, {f\T\i) = 0. 

Let |0) = col(l,0,0), |1) = col(0, 1,0), |2) = col(0,0,l), then we can write the transfer 
matrix in the following form: 

1 

T = 1 • (2) 

1 1 1 

Then the number of states starting with \<7i - ■ ■ cri+Ar_2), and ending with {cTj ■ ■ ■ in 
one limit circle is {ai - ■ ■ ai+N_2\T^^~'^^\aj ■ ■ ■ 

The transfer matrix contains all the information about the time average of space de- 
pendent functions, such as the one point function (cij), and the two point function (aiaj). 
Here (■ ■ ■) means the average in time over a limit circle. Define I = { the complete set 
of states |o-L-Ar+2 • ■ ■ allowed by the boundary condition (H) }, ^ = { complete set of 
states |(7i ■ ■ -(Jiy^i) in general }. The one point function can be evaluated in terms of the 
summation of all states in a limit circle, i.e. 
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where we have introduced the "evaluation" matrix E = J^rn \Vi)'^i{Vi\ signify the fact that 
the sum over rji can be expressed as a matrix independent of the position i. Similarly, the 
two point function is 

1 



Vi,Vj a(^A fSGl 
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If the power of T can be evaluated, then it is easy to get these functions. The power of T 
can be evaluated by the diagonalization of T . Alternatively, the direct multiplication of T 
will show iterative relations which maybe easier to solve and will be used in the following 
section. 



III. N = 2 SANDPILE MODEL 



For the N = 2 model, T" has the general form 



rpr. 
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The iterative relation can be obtained by T""^^ = TT", it is a^+i = 6„ and = 2a„ + hn- 
Therefore both of them satisfy the iterative relation 



Un+l = Un + 2Un-l. 



(6) 



The initial condition of these equations can be calculated from T,T'^,T^ by direct multipli- 
cation. They are ai = , 02 = 1 , as = 1 and 61 = 1 , 62 = 1 , &3 = 3. The general form of 
an and bn are a„ = = | (2"-^^ — (—1)"'^"'^^. Substitute a„ and 6„ into @,(^, the results 
are (for i < j) 



(a,) = 3/2 + (-ir2--i 

= 9/4 + T'^-^ + 3 (-1)"' 2''-^ + (-1)^' 2'^-^ 

Then two point correlation function is (for i < j) 



(7) 
(8) 



These results agree with however the derivation here is much simpler and more general. 
The three point function can be evaluated in a similar way. We got the three point function 
(for i < j < k) 



{a.ajak) = 3 (-1)-^+'^ 2^-''-^ + (-1)"'= 2-''-^ + 3 (-1)"'+'' 2'-^-^ + 9 (-1)"' 2"'-=^ + 

3 2-^-^ + (-1)-'+' T-''-' + (-1)*-^'+'= 2-'+^-^-2 + 27/8. (10) 

The three point correlation function is (for i < j < k) 

Noting that hi = Y.k=i+i'^i^ ^^so get the L dependence of the correlation function of 
height exactly as follows 
lA/,,2, 5L 1 2 1 2 1 8 L 



'"^ 12 36 5L 9-22^' 9L • 22^ 9-24^ 45L ■ 2^^ 6 ■ 22^^ 

Z=l 

(12) 

For large L, (1/L) Eti ((^f) - (^i)') « L, which agrees with what is expected qualitatively 
from the random walk argument in |^. 

IV. iV = 3 SANDPILE 

For the case of = 3, the situation is a little more complicated. The size of the matrix 
T become larger. T is a 4^ by 4^ matrix and the calculation is somewhat more elaborate. 
The process of solving the case of = 3 will show more structure of the transfer matrix 
which will be useful when dealing with the cases of > 3. The elements of T for N = 3 
are {ab\T\cd) = 6bc (1 — (1 — '^3a)(l — ^36)(1 — ^sd)) by the constraint (|l|). 

There is an obvious block structure of T for = 3 model. In order to represent this block 
structure succinctly, we shall adopt the following basis which can be generalized to arbitrary 
A^ later. Let bi be the column vector with component (bi)^ = 6ij. Assign bi ■ ■ ■ bg to \ab), 
bio • • ■ bi2 to |a3), bis ■ ■ ■ bis to |36), and bie to |33), where a,b = 0,1,2. The order of 
blocks is given by the successive sequence of a binary expression if one replaces 3 by 1 and a,b 
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by 0. The relative order between different \ab)s or |a3)s, |3a)s is not important, but we take 
the order as the number sequence in base 3 for convenience. Let A = {\ab)\a, b = 0,1, 2}, 
B = {\a3)\a = 0, 1,2}, C = {\3a)\a = 0, 1,2} ,and D = {|33)} to be the subspace of original 
space spanned by whole vector in the basis. 

The block structure for N — 3 model becomes explicit in T", n > 3, i.e. , the matrix 
of composes of the "saturated' block. A matrix is called "saturated" if all the elements 
are equal. Let E"*^" be the "saturated" matrix of m rows and n columns with all elements 
equal to 1. Then T" can be expressed as 



A9{\ab)) 


Bs{\a3)) 


Cs{\3a)) 


mm 


A^ilab)) 






bn-l^'""' 


bn-l^'""' 


B,{\a3)) 




a„S3x3 


6nS3x3 


bn^'""' 


Cs{\3a)) 




C„E3X3 


a„S3x3 


a„E3xi 


A(|33)) 
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where an,bn,Cn are called block coefficients satisfying the same iterative relation 

Un-^l = tin + 3lin_i + 9w„_2 (14) 

but with different initial conditions. The initial conditions are 

oi = 0, a2 = 1, aa = 4 
6i = 1, 62 = 1, &3 = 4 

C2 = 1, C3 = 1, C4 = 7. 

The block structure helps in calculation. If the trivial E"^^" matrix is ignored, we can use a 
4 by 4 matrix instead of 16 by 16 to represent it. This is the basic idea of 'reduced transfer 
matrix' which will be discussed later. 

The one point function and two point function can be worked out in the same way as 
the N — 2 case. They are 

3"' {u' + Ijf) 



{(Ti) = 2 + 
{<^i<^j) = \ ( 



3-'up-'-^ + 2a;*-^' + 3a;-^' + 4cu-*J + h.c.) + 4, 
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(15) 
(16) 



where a; = — 1 + y— 2- The two point correlation function is (for i < j) 



(cr,a,)e = {a,a,) - = - + (2 - 3 



+ h.c . (17) 



Since the structure of the block submatrix in the transfer matrix is simple, the original 
matrix can be reduced to a simpler 'reduced transfer matrix' to represent the iteration 
relation. We use capital letter to label the block, for example, T^b nieans the submatrix 
which maps the subset A of basis into B. The transfer matrix has the following special 
property which produces the block structure in T" for n > 3. The sum of the elements in 
a row of a submatrix, T^^b, is independent of which row one sums over. Define the sum 
of the elements in any of the row of Tab as Tab- Note that we have managed to reduce 
each submatrix Tab to a number Tab- Using these reduced matrices, one can simplify the 
multiplication between T and T" {n > 3) by inventing a reduced matrix i?" corresponding 
to each T" by ignoring the trivial S matrix. Denote the block coefficient of Tj^ as Rjk- 
Note that i?" is a 4 by 4 matrix for = 3 case. One can easily show that = TR". 

The T for = 3 model is 



T 



10 

3 1 

3 10 

3 1 



(18) 



Multiply Eqn.(|TB|) and the reduced matrix corresponding to Eqn. (|TB|) we can get the 
iterative relation between a„,6„ and c„ easily. The reduced transfer matrix simplifies the 
original 16 by 16 matrix into to 4 by 4 matrix without losing the information of iterative 
relation. By reduced transfer matrix, the iteration between a„, 6„, c„ can be written as 



fln+l = On + 3c„ 

bn+1 = 3ari + bn 

Cn+1 ~ 3cin_i + (Xn 
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It can be shown that an,bn, Cn satisfying Eqn.(|^). To solve this equation, we have to solve 
the polynomial equation = + 3a; + 9. 

For arbitrary A^, the reduced transfer matrix is a 2'-^"^^ by 2'-^"^^ matrix, and they take 
the form 



Af ■ ■ ■ A^^_, 

Af A^ A^ ■ • • A^^_, 



(19) 



where the Af and A^^ are submatrices of the reduce transfer matrix with 2 columns and 
2Af-2 j^Q^g^ (Af)pq = Sip{N6gi + 5q2)- {A'/^)pg = Sip{6q2)- The eigenvalues of reduce transfer 
matrix plays an important role in solving the model. It can be proved that the eigenvalues 



of reduced matrix are roots of the equation = 



1=0 



N-l- 



\ The correlation function 



can be expressed as the combination of these eigenvalues. This work is in progress. 



V. DISCUSSIONS 

We have introduced the concept of transfer matrix into the study of steady properties 
of sandpile models. This is closely related to the Hamiltonian formulation of the usual 
statistical mechanics. Indeed, one may regard the formulas for the correlation function such 
as equation (§) and @) as a "path integral" expressions in a discrete formulation. The 
transfer matrix plays the role of the evolution operator. 

The usefulness of the transfer matrix formulation was illustrated by deriving the one- 
, two- and three-point correlation functions for a deterministic sandpile model with the 
critical slope iV = 2 as well as the one- and two-point functions for the same model with 

= 3. In the latter case, we found that the two-point function decreases exponentially as 
the separation of the two point increases with a correlation length of in the unit of 

lattice spacing. For the N = 2 case, the correlation length is (ln2)~^. It will be interesting 
to see if the correlation length become infinity in the large N limit so that the self-organize 
criticality in the spatial correlation is restored in the limit. 
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